Import

Covid19 Data

Covid19 Japanの公開データ（2020-10-13付）を使います。

　

Area Data

都道府県と地方区分に関するデータは自作データで都道府県コード（pcode）順で各区分を因子化してあります。

　

Data Wrangling

Summarize

オリジナルのデータがどのようになっているかskimrパッケージを用いてサマライズしてみます。元がJSON形式なので、読み込んだ直後は殆どの変量（フィーチャー）が文字型になっているのと欠損値が多いことが分かります。

df %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 91611
Number of columns 23
_______________________
Column type frequency:
character 19
logical 3
numeric 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
patientId 0 1.00 1 8 0 90154 0
dateAnnounced 0 1.00 10 10 0 259 0
gender 10345 0.89 1 1 0 2 0
detectedPrefecture 0 1.00 3 15 0 49 0
patientStatus 87807 0.04 8 23 0 8 0
notes 47561 0.48 1 270 0 41414 1
mhlwPatientNumber 91162 0.00 1 11 0 434 0
prefecturePatientNumber 8446 0.91 5 20 0 83156 0
prefectureSourceURL 60396 0.34 5 224 0 3428 0
residence 18181 0.80 1 38 0 1414 0
sourceURL 637 0.99 1 239 0 7289 0
relatedPatients 82052 0.10 2 259 0 5858 0
knownCluster 89234 0.03 3 88 0 225 0
detectedCityTown 68221 0.26 2 22 0 661 0
cityPrefectureNumber 68328 0.25 1 34 0 23274 2
citySourceURL 79857 0.13 9 317 0 3605 0
deceasedDate 89966 0.02 10 10 0 209 0
deceasedReportedDate 90502 0.01 10 10 0 179 0
deathSourceURL 90598 0.01 14 123 0 621 0

Variable type: logical

skim_variable n_missing complete_rate mean count
confirmedPatient 0 1 0.98 TRU: 90153, FAL: 1458
charterFlightPassenger 91597 0 1.00 TRU: 14
cruisePassengerDisembarked 91600 0 1.00 TRU: 11

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
ageBracket 0 1 33.47 23.21 -1 20 30 50 100 ▃▇▅▂▁

　

Transform

各変量（フィーチャー）を適切な形式に変換し、地域区分でも分析できるように都道府県データと結合します。

x <- df %>% 
  dplyr::mutate(dateAnnounced = lubridate::as_date(dateAnnounced),
                ageBracket = forcats::as_factor(ageBracket),
                gender = forcats::as_factor(gender),
                patientStatus = forcats::as_factor(patientStatus),
                residence = forcats::as_factor(residence),
                cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE)) %>% 
  dplyr::select(dateAnnounced, ageBracket, gender, detectedPrefecture,
                patientStatus, knownCluster, cluster) %>% 
  dplyr::left_join(prefs, by = c("detectedPrefecture" = "pref"))

x
x %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 91611
Number of columns 13
_______________________
Column type frequency:
character 2
Date 1
factor 9
logical 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
detectedPrefecture 0 1.00 3 15 0 49 0
knownCluster 89234 0.03 3 88 0 225 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dateAnnounced 0 1 2020-01-15 2020-10-13 2020-08-07 259

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
ageBracket 0 1.00 FALSE 13 20: 22532, 30: 13970, 40: 11422, -1: 10435
gender 10345 0.89 FALSE 2 M: 45809, F: 35457
patientStatus 87807 0.04 FALSE 8 Dec: 1635, Hos: 1266, Hom: 316, Dis: 283
pcode 1071 0.99 FALSE 47 13: 28431, 27: 11280, 14: 7680, 23: 5667
都道府県 1071 0.99 FALSE 47 東京都: 28431, 大阪府: 11280, 神奈川: 7680, 愛知県: 5667
八地方区分 1071 0.99 FALSE 8 関東地: 47739, 近畿地: 18115, 九州地: 10277, 中部地: 9146
広域圏 6425 0.93 FALSE 8 首都圏: 47946, 近畿圏: 17568, 中部圏: 7760, 九州圏: 7426
通俗的区分 1071 0.99 FALSE 11 関東: 47739, 関西: 17568, 東海: 7430, 九州: 7426
fct_pref 1071 0.99 FALSE 47 Tok: 28431, Osa: 11280, Kan: 7680, Aic: 5667

Variable type: logical

skim_variable n_missing complete_rate mean count
cluster 0 1 0.03 FAL: 89234, TRU: 2377

　
文字型を因子型に変換するだけでも、変量（フィーチャー）内の水準の比率が大まかにつかめるようになります。
例えば、年代別の陽性判定者数は20代が最も多く、続いて、30代、40代と働き盛りの世代に多いことが分かります。また、都道府県別では、東京、大阪、神奈川、愛知と成っていますが、地方区分別では、以外にも関東、大阪の次に九州がきていることが分かります。

　

Visualizing

全体傾向

sec_scale <- 50

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n),
                ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)) %>%
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity",
                      fill = "dark green", alpha = 0.5) + 
    ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
                       colour = "dark green") +
    ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark red") + 
    ggplot2::scale_y_continuous(
      name = "陽性確定数（日別）・7日間移動平均（濃赤折線）",
      sec.axis = ggplot2::sec_axis(~ . * sec_scale,
                                    name = "累積陽性確定数（折線）")
    ) +
    ggplot2::labs(title = paste0("Covid19, Japan（@", Sys.time(), "）"),
                  subtitle = "", x = "日付")
## Warning: Removed 6 row(s) containing missing values (geom_path).

# plotly::ggplotly()

地方区分別

都道府県別では区分が多すぎるので八地方区分別の積上げ棒グラフを描いてみます。サマリで見たように九州地方が以外にも多いことが読み取れます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

折線グラフで表示してみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::scale_colour_brewer(palette = "Set1")

　
さらに分かりやすくするために関東、中部、近畿、九州の四地区を除く他の地区をその他としてまとめてみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::scale_colour_brewer(palette = "Set1") +
    ggplot2::facet_wrap(~ `八地方区分`, ncol = 2)

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `八地方区分`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "八地方区分", values_to = "n") %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_relevel(`八地方区分`,
                                              "北海道地方", "東北地方",
                                              "関東地方", "中部地方",
                                              "近畿地方", "中国地方",
                                              "四国地方", "九州地方")) %>%
  dplyr::group_by(`八地方区分`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `八地方区分`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19（@", Sys.time(), "）"),
                  y = "累計人数")

　
九州が8月頃から急増しているのは、県別に見ると福岡と沖縄での急増が原因と分かります。

x %>% 
  dplyr::filter(`八地方区分` == "九州地方") %>% 
  dplyr::group_by(dateAnnounced, `都道府県`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `都道府県`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "都道府県", values_to = "n") %>% 
  dplyr::mutate(`都道府県` = forcats::fct_relevel(`都道府県`,
                                              "福岡県", "佐賀県", "長崎県",
                                              "熊本県", "大分県", "宮崎県",
                                              "鹿児島県", "沖縄県")) %>%
  dplyr::group_by(`都道府県`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `都道府県`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19, 九州地方（@", Sys.time(), "）"),
                  y = "累計人数")

　

クラスタ比率

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(`八地方区分`, cluster) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::pivot_wider(names_from = cluster, values_from = n) %>% 
  dplyr::mutate(ratio = (`TRUE` / (`TRUE` + `FALSE`) * 100) %>% round(1)) %>% 
  dplyr::rename(`地方` = `八地方区分`,
                `非クラスタ感染者` = `FALSE`, `クラスタ感染者` = `TRUE`,
                `クラスタ比率[%]` = ratio)

Modeling

時系列（TS）分析

日毎の陽性判定者数を時系列データ形式に変換します。判定者数がゼロの日は報告されていないので、最初の陽性判定者が出た日から日日次のカレンダーデータを作成して結合してます。なお、時系列データ形式の周期については日次データなので7日間を周期として設定しておきます（考え方によっては約1ヶ月ごとで4週間=28日を周期にするのもありかと思います）。

ts_x <- x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
  dplyr::select(-date) %>%
  ts(frequency = 1)
  
tsw_x <- x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
  dplyr::select(-date) %>%
  ts(frequency = 7)

tsm_x <- x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
  dplyr::select(-date) %>%
  ts(frequency = 28)

時系列データに変換したものをプロットすると可視化の項でプロットした棒グラフと同じ形のグラフになることが分かります。

ts_x %>% 
  plot()

tsw_x %>% 
  plot()

tsm_x %>% 
  plot()

上記からトレンド（長期的傾向）を除いたグラフ。デフォルト指定なのでlag = 1

ts_x %>% 
  base::diff() %>% 
  plot()

tsw_x %>% 
  base::diff() %>% 
  plot()

tsm_x %>% 
  base::diff() %>% 
  plot()

トレンド、季節変動（周期変動）、非周期変動に分解した場合。frequency = 1では分解できない点に注意。

tsw_x %>% 
  stats::decompose() %>% 
  plot()

tsm_x %>% 
  stats::decompose() %>% 
  plot()

スペクトル密度の推定。何を意味するのかよく分からない。

ts_x %>% 
  spec.pgram()

tsw_x %>% 
  spec.pgram()

tsm_x %>% 
  spec.pgram()

自己回帰によるスペクトラム解析

spectrum(ts_x, method="ar")

spectrum(tsw_x, method="ar")

spectrum(tsm_x, method="ar")

Infer

ARIMA Model

ARIMA（Auto Regressive Integrated Moving Average, 自己回帰和分移動平均）モデルによる予測を行ってみます。予測に必要なパラメータはステップワイズにより自動的に最適なものが選択されます。

ts_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsw_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsm_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()